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ABSTRACT 



We present an efficient algorithm designed for and capable of detecting elongated, thin features such as lines and curves in astronomical 
images, and its application to the automatic detection of gravitational arcs. The algorithm is sufficiently robust to detect such features even 
if their surface brightness is near the pixel noise in the image, yet the amount of spurious detections is low. The algorithm subdivides the 
image into a grid of overlapping cells which are iteratively shifted towards a local centre of brightness in their immediate neighbourhood. It 
then computes the ellipticity for each cell, and combines cells with correlated ellipticities into objects. These are combined to graphs in a next 
step, which are then further processed to determine properties of the detected objects. We demonstrate the operation and the efficiency of the 
algorithm applying it to HST images of galaxy clusters known to contain gravitational arcs. The algorithm completes the analysis of an image 
with 3000 x 3000 pixels in about 4 seconds on an ordinary desktop PC. We discuss further applications, the method's remaining problems and 
possible approaches to their solution. 

Key words, gravitational lensing - methods: data analysis - techniques: image processing 



1. Introduction 

Gravitational arcs are an important diagnostic for the innermost 
mass distribution in galaxy clusters, and thus they are an impor- 
tant indirect diagnostic for a variety of questions in cosmology 
and structure formation. So far, they have been detected almost 
exclusively by visual inspection of cluster images, although al- 
gorithms for their automated search have recently been pro- 
posed (Lenzen et al. 2004; Horesh et al. 2005; Alard 2006). 

There are several good reasons to search for ways to detect 
arcs in an automated fashion. First, the unambiguous definition 
of arc samples calls for an objective and reliable way to de- 
tect arcs and quantify their properties. Second, wide-field sur- 
veys such as the Canada-France-Hawaii Legacy Survey com- 
bine huge data fields with sufficient depth to reach below the 
detection limit for arcs with their typically low surface bright- 
ness. For arc statistics and its potential importance for cos- 
mological studies or investigations of cluster dynamics, objec- 
tively searching for arcs in these wide-field images promises an 
important step forward. However, scanning wide-field data by 
eye covering hundreds of square degrees for arcs with their typ- 
ical widths of < 1" and lengths of < 10" appears as a hopeless 
endeavour. 

Automated arc searches can be conducted where clusters 
have previously been identified, e.g. through their optical ap- 
pearance or X-ray emission, but they can and should profitably 
be extended to blind searches on large areas. Obviously, such 



goals can only be pursued if an algorithm for automated arc 
detection is available. 

The surface brightness of gravitational arcs is typically 
close to the background, which may vary across the image. This 
is one reason why we propose a new algorithm here rather than 
using one of those that were described and implemented earlier. 
For instance, the anisotropic diffusion underlying the algorithm 
by Lenzen et al. bears the risk of creating elongated features 
from the noise which may then be hard to distinguish from real 
arcs. Moreover, we aim at an algorithm which is simple, thus 
presumably robust, and fast enough to be applied to large data 
fields. This rules out better studied techniques for the detection 
of line-like features in images, such as the Hough transform. 
Finally, the algorithm must be capable of distinguishing arti- 
facts such as diffraction spikes or parts thereof from arcs. As 
we shall show below, the algorithm proposed here does indeed 
satisfy these criteria. 

2. Approach 

We begin the description of our algorithm by summarising 
what it is supposed to achieve. Next, we shall explain the nu- 
merical methods used and the reasons for employing them. 



2.1. General description 

We want the algorithm to identify features in astronomical im- 
ages with the following properties: 
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Fig. 1. Small section of an HST/WFPC2 observation of Abell 2390 
showing three arcs. Considerable substructure is visible in the straight 
arc on the left side. 

1 . having a higher mean pixel intensity than the surrounding 
background, including faint features with a mean pixel in- 
tensity well within the domain of the pixel noise where seg- 
mentation by an intensity threshold is impractical; 

2. being much longer than wide, more precisely being elon- 
gated in the sense that it is possible at each point within a 
feature to find a local direction along which the feature's 
intensity changes much less than perpendicular to it; 

3. being extended in the sense that the feature's local curva- 
ture, i.e. its change of orientation per unit length along its 
principal direction, is small. 

According to this definition, a ring with a circumference of 
100 units (e.g. pixels) and a width of 5 units should be detected 
as one feature, the equally narrow sides of an equilateral tri- 
angle of comparable size should be detected as three features 
because of the high curvature at its vertices, and an approxi- 
mately circular object with a length-to-width ratio near unity 
should not be detected at all. 

Bearing strong gravitational lensing in mind, such features 
may be considerably substructured, which argues for choos- 
ing a detection scale limiting the size of detections to features 
exceeding a certain length, while structures on smaller scales 
should be ignored. Another reason for limiting the detection 
scale from below is the contamination of images with numer- 
ous elongated features generated by noise which are typically 
only a few pixels long. 

A straightforward approach to remove small-scale struc- 
tures is to initially convolve the image with a sufficiently large, 
isotropic kernel function, and to apply all further computations 
to the convolved image. However, while the convolution en- 
hances the signal-to-noise ratio per pixel, it tends to blur and 
remove faint and thin features, besides being computationally 
expensive. It is interesting to note, and easy to see in a Fourier 
representation, that in an image convolved with a Gaussian ker- 
nel, fluctuations of decreasing size become increasingly un- 
likely. 



We thus follow a different approach. We first directly com- 
pute the particular image property we are interested in as local 
average within areas corresponding to the detection scale. This 
avoids the blurring effect of an initial Gaussian convolution, 
but retains its positive aspects, in particular the suppression of 
substructures and the enhancement of the signal-to-noise ratio. 
It also significantly decreases the execution time because aver- 
ages are only calculated for a limited set of pixels instead of all 
image pixels, assuming slow variations of the image property 
spatially close to the pixels in the set. Since they are interde- 
pendent, both the particular image property chosen for the al- 
gorithm and the selection of the set of pixels will be detailed 
later. 

Since different local image properties, for example the 
smoothed intensity or its gradient, will typically vary on differ- 
ent scales, a smooth spatial variation on a single scale cannot 
generally be assumed for all local image properties one may be 
interested in. For example, while the intensity gradient will turn 
by 180° when the position is shifted by a few pixels to oppo- 
site sides of a local intensity maximum, the smoothed intensity 
will hardly be affected by the same spatial shift. The specific 
characteristics chosen for detecting features must thus be taken 
into account in the original selection of the pixel set. 

In the following, pixels combined in a set, together with 
their respective neighbourhoods used to compute local aver- 
ages of certain image properties, are called cells. The cells' 
centre coordinates vary and initialise the pixel set prior to each 
step of the algorithm. Even when averaged within a single cell, 
local image properties turn out to be insufficient for identifying 
faint features in presence of noise. Thus, analogous to lines of 
magnetic flux visible in the pattern of iron filings sprinkled on a 
slab of glass, the similarity of properties among groups of cells 
in the presence of a feature is used as the main criterion in our 
algorithm. 

2.2. Ellipticities and cell transport 

The average local pixel intensity of a feature can vary along its 
length, with the only limitation that it must exceed the bright- 
ness of the background. Thus, it is clearly not a good image 
property to use, even if most of the background could be ig- 
nored after a preselection of possible features above a minimal 
brightness threshold. Besides, this is problematic in itself be- 
cause valid features may be fainter than typical variations in 
the image background. 

However, since the local curvature in the direction of a valid 
feature is small by definition, and a fiat background has no pref- 
erential orientation, directions in the local brightness pattern 
are a sound criterion for a detection based on local correlations 
between image pixels. Since the brightness gradient points into 
opposite directions on opposite sides of the feature, spatially 
averaging over it will eliminate the signatures of faint features. 
Using either the gradient or the structure tensor, which is the 
Cartesian product of the gradient with itself, is thus a prob- 
lematic detection method for thin features and reasonably large 
scale sizes. 
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Instead, one can compute the direction from the ellipticity 
given by the second brightness moments. Let I{x) be the inten- 
sity at the position x = (X\,x-i) and q an appropriately chosen 
weight function. Then, the weighted centre-of-brightness in an 
area A is 

j A Xq(I(x))d 2 X 



j; <?(/«) d 2 * 



(1) 



For example, q(I(x)) = I(x) will return the unweighted centre 
of-light. The tensor of second brightness moments has the com 
ponents 

J A (xi - x t )(Xj - Xj)q(I(x)) d 2 x 



Qa = 



j A q(I(x))d 2 x 



with ;, j e {1,2}. The complex ellipticity is 



X = 



Qn-Qi2 + 2iQi2 
Gn + Q22 



(2) 



(3) 



If I(x) has elliptical isophotes in A with an axis ratio of 
r < 1, then \ - (1 ~~ r2 )/(l + f 2 ) exp(2i#), where 1? is the 
orientation of the major axis of the elliptical isophotes relative 
to the xi axis, x is invariant under rotations of n, as it should be 
because they leave ellipses unchanged. 

A vector e pointing into the direction of the major axis can 
be obtained by bisecting the phase angle of x (see Appendix A), 

, ((Xi + lxlXi) for xi >0 d 
a = < and e = — . (4) 

{(X2,\x\-Xi) for Xi <0 \d\ 

The orientation obtained in this way is more stable against 
changes in the scale size than gradient-based methods and has 
a higher signal-to-noise ratio because the second brightness 
moments are typically computed near a feature's centre-of- 
brightness. 

For recognising multiple features in an astronomical image 
without a priori information on their positions, it is obviously 
necessary to process the entire image initially. Our algorithm 
starts with cells evenly distributed and convering the image 
completely. Higher sensitivity can be achieved if neighbouring 
cells overlap, but a large overlap should be avoided to reduce 
the number of necessary operations on the image. 

The first step in measuring the ellipticity is finding the cen- 
tre of brightness as in Eq. (1). For a feature at or near the edge 
of a cell, the centre cannot be accurately determined in a sin- 
gle computation of x. For example, when applied once to a 
smoothed vertical line with a Gaussian brightness profile on an 
already estimated background Iq, 



I(x) = /(xj) = 7 + e 



-x 2 J2a- 2 



(5) 



with a weight function q(I) = I - Iq and an area A extending 
from to 4cr in the x\ direction, Eq. 1 yields 

jci = xie-^dxij |JT e-*? /2 ^djciJ 

4a- ( l~ 4<T \ l 

= -(crV^^I ) crJ^erf(-^) * 0.798cr . (6) 



The result improves if x is computed several times, each time 
shifting the cell's centre to the position x found in the last it- 
eration. In our example, this yields approximate centre posi- 
tions x h0 = 2cr, x u w 0.798o-, x u ~ 0.210o% Jc 1>3 w 0.048cr, 
jci,4 ~ 0.01 lcr and £1,5 ~ 0.002cr after five iterations, where the 
second subscript denotes the iteration number. 
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Fig. 2. Another section from the HST/WFPC2 image used for Fig. 1, 
now centred on the optical emission of Abell 2390, showing cell paths 
leading to the closest local centres of brightness. The faint diagonal 
line crossing the image from the bottom left to the top right is a diffrac- 
tion spike from a nearby foreground star. 



In order to allow meaningful later correlation measure- 
ments between cells, shifting initially neighbouring cells to ex- 
actly the same final position must be avoided, since there the 
ellipticities would then necessarily be equal regardless of a fea- 
ture's real properties. Ideally, the cells should instead be dis- 
tributed equally along the feature's main axis at the end of their 
iterative motion. More precisely, they should end up on or near 
the feature's ridge line, i.e. the line parallel to its local direc- 
tions along the brightness maxima measured perpendicular to 
these directions. 

Carrying out too many iteration steps may create paths con- 
verging on local maxima due to substructures in the feature, 
while too few may leave x far from the ridge line. In the worst 
case, this may result in an orientation e perpendicular to the fea- 
ture for cells in areas characterised by a positive second deriva- 
tive in the intensity profile perpendicular to the ridge line such 
as "valleys" between neighbouring maxima. This is because 
the orientation e tends to be perpendicular to the direction of 
maximum intensity curvature (the secondary brightness mo- 
ments Qij vanish for an intensity linearly dependend on both 
coordinates). In the case of a Gaussian intensity profile, the ef- 
fect would be visible if cells ended up more than one cr away 
from the ridge line. Thus, some compromise has to be achieved. 
Typically, about three to four iterations have the desired effect. 
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Fig. 3. Cell orientations are displayed in the same area shown in Fig. 2. 
Bold white lines mark orientations with c' > c th . 



It is important to note that this method is invariant under a 
linear scaling of the relevant weighted intensities. In the cur- 
rent implementation of the algorithm, q(I) = max(7 - 7, 0) is 
used for the determination of the centre-of-brightness, where 
7 is the mean intensity in A, leaving the method also invariant 
under a constant offset to / and, more importantly, leading to a 
reasonably fast convergence towards the ridge line. 

2.3. Correlation of neighbouring cells 

Once the cells have been shifted to their new centre positions x, 
normalised directional vectors e can be readily calculated using 
Eqs. (2) through (4), where the weight function q(I) - I - 7 is 
used for determining the second brightness moments Qij. Cells 
not initially near a feature and on a mostly flat background 
move randomly by short distances due to the noise, and their 
orientations will also be random, while cells originally located 
near a feature will arrive at different positions along its ridge, 
with very similar local orientations. 

In the following, superscripts will enumerate the cells. In 
one dimension, the cells' relative spatial ordering is obviously 
preserved: if x' Q < xi, then x' n < x J „ after any number n of iter- 
ations of the centre computation. This is not generally true in 
two or more dimensions, but the number of exceptions remains 
negligible for small n. Therefore, the set N of neighbouring 
cells of a cell at the end of its path is assumed to be the same 
as at the beginning of its path. 

Both spatial and directional information should now be 
used to determine whether a cell is likely to be part of a valid 
feature or not. To this end, we study the correlation of a cell i 
with its neighbours j € N. 

The similarity of orientations between cells can be ex- 
pressed by 



The relative coordinates of the cell j in an orthonormal coordi- 
nate frame centred on cell i with its x\ axis parallel to e' are 



A' 



x' and A 



AV 
-AV 



+ A> 
+ A' 



2' 2 



(8) 



2*1 



The coordinate A2 measures the distance cell j from a possi- 
ble feature through the centre of cell i pointing towards e'. For 
convenience, we introduce a measure in the range [0,1], 



1 _ \M 

■\o~* 



for |A 2 | < d\ 
else 



(9) 



where d\ is the initial distance between neighbouring cells. 

Using only the orientation and distance correlations cl and 
c'J, it is already possible to estimate the correlation c l> = c'Jc'^ 
between the cells 2 and j reasonably well. Including the entire 
neighbourhood, this becomes 



(10) 



where \N\ denotes the number of cells in the neighbourhood N. 

To reduce the effect of possible errors due to the assumption 
of an identical initial and final neighbourhood, and to lower 
the weight of closely neighbouring cells, N can be extended 
and another factor depending only on |A| can be introduced 
which is small for |A| = 0, tends to zero for |A| » d\, and 
has a maximum in (0, d{\. For a sufficiently large initial N, the 
modified correlation measure 



JJJj 



(ID 



cj = \e'-e j \ 



(7) 



becomes independent of the initial neighbourhood and allows 
suppressing the contribution to the correlation of cells very 
close to cell 2. Since it slightly increases the execution time, 
introduces an additional degree of freedom and leaves the suit- 
able form for to be determined, this distance weighting was 
ignored here. 

Values c' are computed from (10) for each cell, restricting 
the neighbourhood to the eight surrounding cells. Cells with 
c' < c t h are excluded from the further processing. For many 
applications, a reasonable value for the threshold is c± = 0.5. 
In a next step, those cells above the threshold which are likely 
part of a feature are grouped into objects, using a very similar 
method as described above. 

2.4. Object generation 

Each remaining cell i is now individually compared with all 
other remaining cells j in its neighbourhood, based on the 
known correlation coefficients c 1 , and ci. Neighbourhoods 
are extended essentially with a slightly modified version of 
Bresenham's line-drawing algorithm applied to cells instead of 
pixels. It generates an appropriately oriented region, for exam- 
ple with a length of nine and a width of three cells, surround- 
ing the cell i in its initial position. To avoid cells not following 
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a possible curvature in the feature, cells with a large separa- 
tion Ai can additionally be handicapped, e.g. by introducing a 
factor 

0.2 



c'J = 0.8 + 



1 + 

1 T 200 \ d, ) 



into the correlation coefficient 



r<i] - 'J >J r K 
~ d x c 



(12) 



(13) 



between any two cells i and j. If C' exceeds a threshold C±, 
cell j is added to the object which cell i is already belonging to. 
If cell i is not part of an object yet, an object is created including 
both cells. 

Cells may belong to more than one object, which is rea- 
sonable because objects cannot be properly distinguished be- 
fore they are completely defined. When the previous step is 
completed, any two objects with at least one common cell are 
combined into one. Since any averaging between neighbour- 
ing cells was avoided which may not be part of the underly- 
ing feature, the value of C± may exceed c t h- Tests showed that 
C t h = 0.7 is a reasonable choice. 

Cells located at the exact same position, e.g. at the centre of 
a radially symmetric feature, will be grouped into one object. 
Since only features exceeding a certain length are valid, such 
objects can later be removed from any following analysis. Also, 
objects consisting of only very few cells are likely to be the 
result of the random behaviour of cells far from any feature, and 
can similarly be invalidated. For testing, we set the minimum 
number of cells in an object to four. 

2.5. Object graphs 

The above procedure may result in objects covering several fea- 
tures. For disconnecting them, the topographical structure of 
the feature, or of spatially connected features underlying them, 
must be taken into account. Defining objects by concatenation 
of cells is not practical for this purpose. Therefore, it is useful 
to represent their ridge lines by undirected geometric graphs. 

2.5.1. Initial graph generation 

For each object, graphs are created in the following three steps: 

1. nodes are defined from spatially close cells, whose aver- 
aged position and orientation are assigned to the nodes; 

2. if another node is found within a limiting angle of ±n/4 of 
a node's direction, it is connected with the node; and 

3. separate connected subsets are connected into one con- 
nected graph. 

For finding the cells to be combined into nodes, each cell is 
initially marked as unprocessed. The algorithm loops over all 
unprocessed cells in the object and again calculates the absolute 
cosine between the cells' directions ci, the position A of cell j 
in the oriented coordinate system of cell i, and a measure of the 
perpendicular distance 



1-3^ for |A 2 |<2c/i 
else 




Fig. 4. Graphs overlaying the final detections. The orientation at- 
tributed to each node is marked by a black and white line. 



for every other unprocessed cell j in the object with a spatial 
distance |A| < d\, where d\ is again the initial distance of the 
cells i and j. 

Using a third correlation threshold CL, the mean positions 
and directional vectors of all cells j with c'J = cjc'* > CL are 
determined, the cells j and ; are marked as processed, and a new 
node with these properties is created. In c''J, the perpendicular 
distance is weighted less than in c'J before to include cells of 
high orientational correlation in a greater spatial range into one 
node, thus preventing the construction of spatially close paral- 
lel nodes which might later cause loops in the graph. 

The algorithm continues by searching the next unprocessed 
cell. It is important to note that the mean of the directional 
vectors does not return an average orientation, which must in- 
stead be calculated as the mean over the normalised elliptici- 
ties, e.g. using 

= ( e i _ el,2e\e 2 ) , (15) 

which can then be used to derive directional vectors. In this 
way, all nodes in the graph are created as cell averages. As 
with C t h, the threshold CL quantifies the relation between two 
instead of more cells, and should thus exceed c t h- Tests showed 
that CL = C t h = 0.7 avoids including substructures but allows 
features to overlap. 

To connect the nodes, the relative distance A from a node i 
to every other node j is calculated using Eq. (8). From this, a 
modified distance 



-V 



Af + 16A 2 2 



(16) 



(14) 



is computed for each node with , i.e. within an angle of n/4. 
The nodes with the lowest r e among all nodes are connected to 
node ;, if they exist. 

The angle tt/4 limits the curvature to less than n/4d\. The 
metric r e assigns equal distances from cell i to all points on an 
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ellipse with an axis ratio of four and a major axis parallel to 
the e\ axis centered on i. Consequently, nodes with a small per- 
pendicular distance are connected preferentially. This reduces 
the influence due to nodes of features crossing within the same 
object. 

The implementation of this method can lead to a discon- 
nected graph, with individual features separated by indepen- 
dent structures in the object. A connected graph can be cre- 
ated by repeatedly connecting the subgraph containing the first 
node in the object to the closest subgraph separated from it. 
Specifically, an additional variable v is assigned to each node, 
which is initially set to zero for all nodes. In a loop, the value 
one is recursively propagated to all connected nodes, starting 
in node zero. In the next step, the distance r e is measured be- 
tween all nodes i connected with node zero which thus have 
v' = 1, and all disconnected nodes j with v' = 0. The pair of 
nodes ; and j with the minimal distance is then connected, all v 
are re-initialised with zero, and the process is repeated until no 
disconnected nodes are left. 



2.5.2. Generation of subgraphs from graphs 



2.5.3. Graph concatenation 

Our algorithm's final step is the concatenation of those new 
graphs which may represent a single feature. To this end, a cir- 
cle is fit to each graph, and the nearest of any other graphs 
falling on that circle is added to it, provided it is closer than 
three times the angle spanned by the two lines from the centre 
of the circle to the first and last nodes in the original graph. 

Then, a new circle is fit to restart the procedure. Since a 
true least-squares fit of a circle to a set of points is impossible 
analytically, and finding the solution numerically is computa- 
tionally expensive, the fit is simplified by assuming that one of 
the first and the last two nodes of the graph fall exactly on the 
circle. The radius can then be fit analytically, and the best fit 
in a least-square sense is chosen from the four solutions. The 
deviation of this method from an unconstrained fit is negligible 
in most cases, although the method is feasible only because of 
the already simple form of the given graphs. 

3. Potential problems 

A number of unresolved problems remain, as well as problems 
with known solutions. Both are considered in the following. 



The graphs resulting from the preceding procedure are not ide- 
ally suited for post-processing because they may belong to mul- 
tiple features, or features contaminated by spurious structures. 
Therefore, all their subgraphs which possibly represent either 
part of or a complete valid feature are generated using the spe- 
cific structure of the graphs representing valid features: each 
node must be connected to either one or two other nodes, and 
the relative angle between two edges at a node must fall within 
7T + a, with < a < tt/4 depending on the maximum allowed 
curvature. That is, graphs satisfying this definition are either 
arc-like, with two end nodes having only one connection and 
an arbitrary number of nodes with two edges in between, or 
ring-like, having no end nodes and containing only nodes met 
by two edges. 

To find the first variant, all edges connected to a node ; 
are checked for edges joining on the opposite side of the node 
within the specified angle range. If there is no such opposite 
edge, node i is a starting or ending point, and the algorithm re- 
cursively goes through all nodes connected through this and the 
following edges which are within the required angle interval 
with their preceding edges, and lead to nodes not already be- 
longing to the current path. All nodes j in this recursive scheme 
for which no successive edge is found within the angle interval 
are end nodes. A new graph is created starting with the path 
from ; to j if an equal graph, consisting either of the same path 
from i to j or its exact reverse, does not exist yet. 

To find the second variant, one can start the recursion for all 
edges in the original graph which are not yet part of any other 
graph. The resulting graphs will miss one edge and thus have 
a starting and an ending point, which makes further computa- 
tions slightly easier, but will equally represent any underlying 
feature. We tested our algorithm setting a - n/5. 



3.1. Ellipticity bias 

Given an odd scale size do, the straightforward method of mea- 
suring ellipticities for each cell computes the sum 



Qij = 



XxeA(Xi ~ Xj)(Xj - Xj)q{I(x)) 



(17) 



inside the rectangle 

A = [x e Z 2 : xf n < x x < x™ x A xf n < x 2 < x m ™ 



with x n 



<h -I 

X\ Z ,X2 



l 2 

do - 1 



2 2 
and x nmx = (xf n +do,x™ in + d ) , 



(18) 



where x are pixel coordinates within the cell. However, since 
noise fluctuations in the four corners of A have a stronger influ- 
ence on the secondary brightness moments than similar fluctu- 
ations along the edges or in the interior of the cell, the resulting 
ellipticities will preferentially be aligned with the diagonal. 

This significant ellipticity bias can be avoided using 
position-dependent weights q'(I(x),x) = q(I(x)) ■ A p (x) and 
modified coordinate components x[(x) and x' 2 (x) in (17), e.g. 

ZxgaC^ ~ *«)(*') - Xj)q'(I(x)),x) 

Qtj = v ^7-—- , (19) 

where A p (x) is the area of overlap between a circle around x 
with radius y and the horizontally aligned square of one pixel 
side length centred on and typically associated with a single 
pixel at position x. The centre of A p (x) is x'(x), and the total 
integration area A is the same as above. Using only discrete 
centre coordinates x, these values can be pre-calculated analyt- 
ically during the initialisation of the algorithm (see Appendix 
Q. 
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3.2. Influence of point sources and galaxies 

Several problems are mostly associated with point sources in 
the astronomical context, although they can also appear in other 
applications. 

For a single bright point source, one problem already solved 
by setting a minimal object length is the creation of objects 
from multiple cells converging to the exact same location. An 
efficient method to approximate the length of an object repre- 
sented by a set of cells is to first find the cell j farthest from an 
arbitrary cell i, and then the cell k farthest from cell j, using the 
distance of cell k from cell j as a length measurement. This is 
not necessarily the largest pairwise distance in the object, but it 
can only be smaller by a factor of V3 (see Appendix B). 

A closely related problem for very bright point sources is 
caused by brightness slopes which may be much more extended 
than the scale size do- Since a fixed number of iterations is used 
to determine the nearest local centre-of -brightness, cells do not 
necessarily complete the path to this point and measure the sec- 
ondary brightness moments on the slope. They will then have 
orientations parallel to the slope's gradient and can, depend- 
ing on the original cell position, generate spurious detections. 
This can be avoided by setting limits for the total spatial dis- 
placement of a cell and discarding cells outside from further 
computations. Using q(I) = max(7 - /, 0) and assuming a con- 
stant slope, the total spatial displacement is the number of it- 
erations times do/3, and less for local brightness maxima with 
their negative curvature in the brightness profile. We set the 
limit to 8<io/ 1 1 for testing the algorithm. A welcome side effect 
of this approach is the preservation of the valid cells' original 
neighbourhoods. 

As a consequence of averaging over an area characterised 
by the scale size, bright pixels in the vicinity of a valid feature 
can outshine it and prevent its detection. Since point-source im- 
ages must be closer to a feature than about half the scale size do 
to pull cells away from it during the centre-of -brightness deter- 
mination or to influence the ellipticity measurement, choosing 
a smaller value for do can avoid this problem at the cost of a 
lower signal-to-noise ratio. 

Clustered point sources can create false detections by im- 
itating substructures of a valid feature. Since the distinction 
of substructure from point-source clustering is impossible af- 
ter averaging over do, which is one of the basic ideas behind 
the algorithm, post-processing of the brightness profiles under- 
lying each object graph is possibly the only way of solving this 
type of problem. Even a single point source enhances the likeli- 
hood for a spurious detection by pulling multiple cells towards 
a single location, thereby compromising the simple cell-count 
filter. 

Adding a single cell j with sufficient correlations c 1 and O' 
further away than the minimal object length and in the neigh- 
bourhood of any cell ; among the already assembled cells re- 
sults in a detection. For farther cells with high spatial correla- 
tion, i.e. placed along the assembled cells' mean direction, the 
orientational correlation will also be increased if they are af- 
fected by the point-source image, causing them to point at the 
source. 



The situation is similar for two point sources, which can 
create false detections if their distance in the image exceeds 
the minimal object length and if a cell located at the centre of 
each point source image can "see" the image of the other within 
its averaging area. Of course, a larger point-spread function in- 
creases the influence of the points sources. 

Searching for arcs, extended diffraction spikes and bloom- 
ing satisfy all criteria for valid features and will thus cause spu- 
rious detections. One approach to avoiding them is to deter- 
mine the location of diffraction spikes in an independent step 
by searching for bright point sources, then masking bloom- 
ing by its maximal brightness and approximating the extended 
point-spread function. If an arc candidate is found which spa- 
tially coincides with a diffraction spike, it can then be flagged 
or invalidated. 

Searching for prominent arcs, the scale size may be too 
large for detecting galaxies, but if the scale size is small 
enough, they will also cause spurious detections if they ful- 
fil the criteria for valid features. They can then be eliminated 
based on their lower length and length-to-width ratio of their 
isophotes. Another possible way to invalidate them is to com- 
pute their radial brightness profiles, which will generally be 
shallower for gravitationally lensed objects. If small arcs are 
included into the search, however, it may be impossible to re- 
liably remove them unless observations in multiple frequency 
bands allow a colour discrimination. 

3.3. Chip and image boundaries 

When image data from CCDs of different resolution and sen- 
sitivity are combined into one image, e.g. for the HST/WFPC2 
instrument, or if the image examined is a mosaic, boundaries 
across the image between regions of different noise level or 
mean background intensity may occur. A noticeable change in 
the mean background intensity shifts cells by about do/2 to- 
wards the brighter region, significantly increasing the cell den- 
sity in this narrow space and possibly creating an orientation 
bias perpendicular to the boundary. Since both directional cor- 
relation and a low spatial distance perpendicular to each cell's 
direction are required for the creation of objects, this case will 
not necessarily result in spurious detections, although they will 
become more likely due to the increase in cell density. 

A change in the noise level and similar mean background 
intensities in both regions typically only introduces an ellip- 
ticity bias in the direction of the boundary, thereby enhanc- 
ing the directional correlation. Using a weight function q(I) = 
max(7 - /, 0) for the determination of the centre-of-brightness 
which ignores pixels with intensities below / additionally leads 
to a cell shift into the region of higher noise and an increase in 
cell density. For these reasons, regions of different noise in a 
single image may cause spurious detections. Setting q(I) = 
for pixels outside the image can be seen as a boundary between 
regions of equal background but with a noise change equal to 
the noise in the image. In the current implementation of our al- 
gorithm, false detections at the image boundaries are avoided 
by simply invalidating cells straddling them. Another approach 
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could generate noise outside the image, with the drawback of 
decreasing the signal-to-noise ratio for the affected cells. 

3.4. Noise and background gradient 

We measure the sensitivity of the algorithm using a model im- 
age of 1 000 x 1000 pixels with intensities drawn from a Poisson 
distribution. The signal is a smoothed line segment of a quar- 
ter circle with 500 pixels radius, randomly oriented and shifted 
by at most 50 pixels from the image centre. Intensities without 
noise are given by I(x) — Io + 7 S exp {^-d 2 (x)j2cr 2 ^, where To is 
the background intensity, / s is the peak signal intensity along 
the ridge of the line segment, d(x) is the distance of the pixel 
at x from the line segment and cr is a smoothing scale. The 
actual intensities in the image are Poisson distributed random 
numbers with a mean of I(x) and standard deviation of y/I(x). 

Although the resulting model signal is not identical to a 
line segment with intensity y/2ncrl s above the background 
smoothed with a Gaussian kernel, it is a good approximation, 
easily constructed for comparison with other algorithms and 
equally close to the typical valid feature in astronomical appli- 
cations as a line segment smoothed with a Gaussian. 

For each set of parameters, we note both the number of 
valid detections N^m, zero or one, and the number of spurious 
detections A^ spui i ous , where a detection was counted as valid if all 
nodes in the detected graph were inside a maximum distance 
of one cr from the arc and the graph was at least 196 pixels 
long, one fourth of the line segment's length. If more than one 
detection fulfilled this criterion, only the first was counted. 

For a constant background of Io = 100, a central signal in- 
tensity of / s = 6 and cr = 10, the mean number of spurious de- 
tections in ten runs for various scaling sizes is given in Table 1, 
showing the detection of short, random structures in the back- 
ground. For the range of scale sizes shown, the number of valid 
detections was zero. 



Table 1. Dependence of spurious detections on the scale size. 
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0.3 


0.1 






For the same values Iq and 7 S , the average detection ratio 
A^ V aiid/(l + A^ S purious) found in ten runs is plotted as a function of 
the scale size do for cr = 6, 8, 10, 12 in Fig. 5. For scale sizes 
below 2cr, the cell's areas are dominated by noise even close to 
the feature, preventing its detection, and for scale sizes above 
10cr, the feature cannot attract enough cells to sufficiently in- 
crease the first correlation value c'. While there are detections 
on the line segment for large scale sizes, these are too short to 
be valid. With increasing cr, the integrated signal intensity in- 
creases proportionally, improving the likelihood for detection, 
as evidenced by the larger range of scale sizes with Alalia = 1 ■ 

For Io = 100, cr — 10 and the scale size do - 45 
best adapted to this cr, the average detection ratio iV V aiid/(l + 
A^spmious) in 50 runs for different signal intensities 7 S is given in 
Table 2. 




Fig. 5. The normalised number of valid detections, N V aiid/(l +-A spu rious), 
is displayed as a function of the scale size d for different smoothing 
scales cr. 

Table 2. Detection ratio for I Q = 100, cr = 10 and d = 45 depending 
on the central signal intensity Z 8 . 
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5.0 
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"valid 
1 +^spurk>us 





0.03 


0.12 


0.38 


0.42 


0.63 


0.89 


1 



In both tests, a noisy but on average constant background 
was assumed. If the background has a slope steeper than the 
intensity slope up to the ridge line for a background-subtracted 
feature, cells cannot find this feature's location during the de- 
termination of their centre-of-brightness, prohibiting detection. 
In astronomical observations, bright foreground objects can in- 
duce steep intensity slopes which are generally not well de- 
scribed by a first-order approximation of the local background, 
thus offering a challenge to possible background-subtraction 
methods. 

4. Implementation specifics 

During the implementation of the above ideas, several arbitrary 
choices were made, some of them with a bearing on the results 
below in 5, others in an effort to make the code more efficient, 
and some to keep the source code manageable. 

4.1. Initial cell placement 

The most important free parameter, and ideally the only one 
to be changed for different image data, is the scale size do- 
Based on it, cells are initially placed on a regular rectangular 
grid spaced by d\ = do/2. Determining the closest centre-of- 
brightness is done by applying a discretised version of Eq. (1) 
to find x in a rectangular area A as defined in (18) for each of 
three iterations. Using a rectangular area and thereby allowing 
farther spatial shifts in the diagonal direction instead of a cir- 
cular disk for this part of the algorithm compensates for the 
initially larger diagonal distance of the cells. 

4.2. Classes and local indices 

Since the source code was written in C++, it was convenient 
to use structures and classes to represent cells, objects and 
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graphs, providing easy access to data members and hiding most 
of the necessary memory management and initialisations from 
the primary functions. Since the number of cells remains con- 
stant while the algorithm proceeds, it was sufficient to imple- 
ment them as one-dimensional arrays of cell structures, where 
the original position and neighbourhood can be readily derived 
from each cell's index. Graphs are described by classes for 
nodes, single graphs, and the list containing all graphs. To in- 
crease the efficiency of the graph concatenation, whose com- 
plexity is 0(n 2 ) in the number of graphs n, the graph list can 
be used to subdivide the image into rectangular regions, each 
of which maintains a list of all graphs with nodes in it. These 
lists can be initialised with linear complexity and are, com- 
bined with a method to mark already processed graphs, used to 
enumerate only those graphs in relevant areas. However, while 
significantly decreasing the execution time, this does not re- 
duce the concatenation's squared dependence on the number of 
graphs. 

4.3. Performance 

Applying the algorithm to a FITS image of 3021 x 3021 pix- 
els with single (floating point) precision using a scale size do 
of 27 pixels and a very low minimal length of 22 pixels on 
a personal computer with one 2.80 GHz processor and 1 GB 
RAM took approximately 4.4 s total execution time of the cal- 
culation thread, and resulted in 137 detected objects prior to 
applying any later filters. Apart from setting a minimal object 
length and cell count, no further filters are included yet. Most of 
that time, namely ~ 3.9 s, was needed for the determination of 
the centres-of-brightness and the orientations. The graph con- 
catenation needed ~ 0.1 s. 

Changing the scale size to 13 pixels and reducing the mini- 
mal object length to 10 pixels resulted in 1033 mostly spurious 
detections of small features, and increased the execution time 
to ~ 1 1.8 s, of which the determination of the centres of bright- 
ness and orientations needed ~ 5.6 s and the graph concatena- 
tion another ~ 5.6 s. 

Applying the algorithm with the latter parameters to an- 
other single-precision FITS image with 8500 x 8300 pixels and 
less small-scale structure gave 395 detections, many of which 
are diffraction spikes, and took ~ 161.2 s in total, 99.9 s of 
which were used for the centre-of-brightness and orientation 
computations, and ~ 3.5 s seconds for the graph concatenation. 
The remaining ~ 57.7 s were used for transferring data between 
the threads and initialisation. Considering that the area is only 
~ 8.5 times larger, the increase by a factor of ~ 25.6 in ex- 
ecution time for the algorithm's first step is mostly caused by 
extensive memory swapping, as is the steep increase in the time 
needed for the data transfer and initialisations. Apart from con- 
ventional optimisation or using more RAM, this can easily be 
solved by applying the algorithm several times on distinct over- 
lapping subregions in the image. 

5. Results 

We developed our algorithm primarily using an HST/WFPC2 
observation of Abell 2390 with an exposure time of 2100s in 



the F814W filter as the base image, which includes multiple 
arcs. We also tested the algorithm against other images to con- 
firm its generalisability. 




Fig. 6. This section of the image was used to "train" the algorithm. It 
is shown here without additional markers. 




Fig. 7. The structured arc in the same image is shown with graphs 
overlaying the detections. The scale size is d = 27 here and in all of 
the other HST images shown. 



6. Future work 

Though it is one aim of the presented algorithm to introduce as 
few as possible arbitrary methods and variables, several exist, 
and they should be optimised. This refers to the free parameters 
already introduced in the text, but also general methods. For ex- 
ample, one possibility would be to change the initially rectan- 



10 



Gregor Seidel & Matthias Bartelmann: Arcfinder 



m 



% - ' 



» v ' 




» ■ • 



m - 



Fig. 8. Another WFPC2 observation of Abell 2218, taken with a 
702 nm filter and an exposure time of 1900 s. 




Fig. 9. This image shows detections marked by black and white lines 
obtained with a scale size of do = 25. Prominent spurious detections 
can be seen near the horizontal border between the WFPC2 CCDs. 



gular grid to a hexagonal one to provide an equally spaced set 
of immediate neighbours and possibly a reduction in necessary 
overlap. 

Point sources can, as detailed in 3.2, significantly increase 
the likelihood of spurious detections, which are probably best 
removed by a post-processing of those areas in the image un- 
derlying detected graphs. Since these areas represent only a 
fraction of the complete image, more elaborate filtering tech- 
niques are possible, but as of now, a filter capable of distin- 
guishing arc substructure from clustered point sources was not 
implemented yet. Stellar photometry in crowded fields already 
deals with similar issues, and trying to model a feature as a su- 
perposition of point-spread functions, similar to the approach 
of the daophot program (Stetson (1987)), could be a sound 
if computationally expensive approach. Using the large elon- 
gation of valid features, the width of low-intensity isophotes 
could also be used to discriminate features induced by only 
one or two point sources, although a slope in the detection's 
background increases the minimal isophote intensity, possibly 
rendering this method useless and necessitating careful back- 
ground subtraction. 

Diffraction spikes are another problem mentioned in 3.2, 
and a method to estimate their positions without specific a pri- 
ori information on the detector is already implemented, but 
must be combined with the our algorithm to remove detections 
coinciding with them. 



We presented a first quantification of the algorithm's sen- 
sitivity and parameter dependence in 3.4, but the simple way 
of modelling a valid feature and the assumption of a constant 
background intensity make it difficult or impossible to extrap- 
olate to the applicability of the algorithm on real datasets from 
this. Also, the observational data presented above was less than 
comprehensive, since comparatively high-quality images were 
used. For these reasons, the algorithm will be applied to simu- 
lated strong-lensing images with varying seeing conditions and 
detector properties as well as to data available from current sur- 
veys, e.g. the HST-based COSMOS or ground-based observa- 
tions which, given their wider area, are particularly interesting. 

As of now, the algorithm uses only images taken in one 
spectral filter. However, where sufficiently deep and resolved 
observations are available in several spectral bands, this infor- 
mation could be used both to remove spurious detections sepa- 
rable by the varying spectra of their components and to increase 
the algorithm's sensitivity by using the additional parameter 
space for correlation measurements. 



7. Summary 

We have presented an algorithm for the fast and automatic 
detection of arc-like features in astronomical data. The algo- 
rithm proceeds in three major steps. First, the image pixels are 
grouped into square cells which are moved in a fixed number 
of iteration steps towards their centre of brightness. They thus 
find a final position on or near a local intensity peak. Second, 
the quadrupole moment of the light distribution is measured in 
each cell, where it defines an elongation and a direction. Third, 
the neighbourhood of each cell is searched for such cells whose 
intensity distribution points into a similar direction. If so, the 
cells are connected to form an object, and the procedure is con- 
tinued until all cells have been processed. The objects found in 
this way can then be automatically classified. 

Having arc-like features in mind which may be just above 
the noise limit, we avoid smoothing in the algorithm which may 
give rise to subtle biases. Smoothing may make arcs disappear, 
but it may also create spurious arc-like features by anisotropic 
stretching of positive noise fluctuations. The algorithm also 
avoids object detection and uses only the local intensity distri- 
bution in the cells. These design criteria enable an implemen- 
tation which operates very fast. For example, it was possible to 
scan an image of 3021 x 3021 pixels in less than five seconds 
on an ordinary desktop PC. 

We used HST images to illustrate the steps of the algorithm 
and its efficiency. Remaining potential problems concern the 
presence of bright point sources and extended galaxies in the 
images, boundaries of images and CCD chips in image mo- 
saics, and residual gradients in the noise or the image back- 
grounds, and we outline how they can be overcome. 

Apart from several parameters which control how cells may 
be connected with neighbouring cells to form objects, the sin- 
gle main parameter is the number of pixels to be grouped into 
a cell. This can be suitably chosen depending on the resolution 
and the quality of the image, and the algorithm is fast enough 
to allow calibration runs. 
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Fig. A.l. An elliptical isophote and the geometric construction of the 
orientation &. 



The main application which we have in mind for the algo- 
rithm is blindly scanning wide-field images for arcs in order to 
construct unbiased arc samples from large data sets. We shall 
study in a forthcoming paper how the detection efficiency and 
reliability will depend on noise, seeing, and the density of fore- 
ground objects. 
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Appendix A: Obtaining an orientation from the 
complex ellipticity^ 

In C, the real axis and the complex ellipticity x — {77? exp(2i#) 
span 2i?, twice the angle 1? in image space between the xj-axis 
and the main axis of an underlying ellipse. Writing Of 1,^2) = 
(Rix), 3 Of)), a vector parallel to the ellipse's main axis can be 
obtained as d\ = Of 1 +\x\,Xt), one diagonal of the rhombus with 
sidelength I = \x\ seen on the right side of Fig. A. 1 . The angle 
between di and the xi-axis is §. Since \di\ — » for2# — > n, this 
method cannot be used for the bisection of angles 2§ close to n. 
Of course, one can also bisect 2ft using di = Of 2, \x\ ~ Xi)> i- e - 
the vector Of 1 — \xVxi) rotated clockwise by |, where Thales' 
theorem using the semicircle in Fig. A.l can be used to show 
that the angle between (x\ - \xVx2) and d\ must be |. Contrary 
to \d\, \d 2 \ -> for 2§ -> 0. Using 




Of 1 + WlX2) for x\ ^ 

(xz,\x\ -x0 for X\ < 



as in (4), \d\ > V2U-I in all cases, making a reliable computation 
of the orientation e = A possible. 



/ 1 
/ / J 

/ / r/ 
/ 1 / 


A 

/ 

/ 

/ 


c^T ^ 

V \ 


/ * / ' 
1 1 / / 

/ '/< 

1 ^ / / 




\ ! c \ 




/ / / X 




f or V 




V 


D 


A l 




\ x \ x 
\ x \ n 
\ x \ x 

\ x V- 

\ x \ ^ 
\ x \ N 


d 


/ 1 /(U 




\ s \ 


s 

\ 

\ \ 

\ s 


\ 1 /I 


\\ 




/\ / 1 
/ / N 

\ 1 



Fig. B.l. The length determination method used in 3.2 requires that 
the centre coordinates of all cells in the object must be in the inter- 
section of the circular disks J{ and 3. It is shown that V3r is greater 
or equal to the maximum pairwise distance of the cells in the object, 
which must be lesser or equal then d if r is between \AB\ (left dashed 
lines) and V2|A6| (right dashed lines). For r > ^fl\AB\ the maximum 
pairwise distance is limited by 2|A6|. 

Appendix B: V3 relation of the object length and 
the maximum pairwise distance 

In 3.2, the following method is used to determine an approx- 
imate object length: First, an arbitrary cell A € Wl, where M 
is the set of all cells in the object, is choosen. Then, the cell 
B € M farthest from A is found. Last, the cell C € M far- 
thest from B is determined where r = \BC\ is used as the object 
length. It can be shown that V3r is greater or equal to the max- 
imum pairwise distance for all cells in the object: 
Putting the method in more formal terms, 

A, B e M with \AB\>\AX\ MXeM and (B.l) 

B, C e M with \BC\ > \BX\ VIeM (B.2) 
=> \BC\ > \AB\ and \BC\ < \BA\ + \AC\ < 2\AB\. (B.3) 

From (B.l) and (B.2), all points in M are inside the intersection 
of the circular disk J{ around A with radius | AB\ and the circular 
disk !B around B with radius r = \BC\. The point where the 
line between the intersections Do and D\ of the boundaries dJi. 
and dS crosses the straight line AB shall be D and the distance 
between Do and D\ shall be d, as illustrated in Fig. B.l. Using 
(B.3), the problem can be separated into three cases: 

1. r = \AB\: J[ and S have equal radius, therefore the 
problem is symmetric and \BD\ = ^. From this, 
(|) 2 = r 2 - |BD| 2 = r 2 - (§) 2 => d = V3r. Since 
j > \BD\ = \AD\ and the circle with radius | centred on 
D therefore contains the complete intersection of J{ and 
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Fig. C.l. Overlap of pixel areas with the area of a circular disk. The 
coordinates (x, y) are chosen so that the disk's centre is at their origin 
and the relevant pixel appears in the upper quadrant enclosed by the 
diagonal lines. To better illustrate the darker areas relevant to the in- 
tegrals (C.l) to (C.4), an even diameter of four pixels was chosen; in 
the algorithm, an odd diameter equal to the scale size d is used. 

S, the maximum pairwise distance must be lesser or equal 
than d and V3r must be greater or equal than the maximum 
pairwise distance. 

2. \AB\ <r< y/2\AB\: a = IDqAB and B = lABD = lBD A 
as displayed in Fig. B.l. For r > \AB\ follows a > | 
and therefore B < |. Using this, cos/3 > \ and 



(!) 



\BD? 



r 2 cos 2 jS < \r 2 => d < V3r. For 



example using \BD\ = \AB\{\ - cos or) and \ = \AB\ sin or 
it is easily shown that j > \BD\ > \AD\ and d must again 
be greater or equal to the maximum pairwise distance. 
Consequently, V3r is greater then the maximum pairwise 
distance. 

3. '\f2\AB\ < r < 2\AB\: the intersection between the circular 
disks J[ and & is equal to J[, therefore the maximum pair- 
wise distance must be lesser or equal to 2\AB\ < V2r < 
V3r. 



for computing parts of the area A p , and 
xdxdy = ^ x(V/? 2 - x 2 -y )dx 



J-(R 2 -x 2 f' 2 -^ 



as well as 
JJ" ydxdy 



n 

p 1 



ydxdy 



(C2) 



(C3) 



(C4) 



(/? z _ x - _ y - Q )dx = -x{R 2 - y 2 Q ) - -x 
2 2 6 Ao 

(C.5) 

for finding the centre x' . The remaining area integrals have 
fixed limits of the form JJ A dxdy = 1 f 1 dxdy and can be 
included in a straightforward manner. 
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Appendix C: Determination of the pixel/disk 
overlap area A p (x) and centre x'(x) 

In 3.1 the intersection A p (x) of a circular disk with radius y 
around x and the square of one pixel sidelength representing a 
pixel at position x must be computed as well as the centre x'(x) 
of this intersection. With a circle radius of R, pixel coordinates 
(x,y) chosen such that the circle's centre is at (0,0), y > 
and \x\ < \y\, and A being the area inside the pixel square, the 
problem mostly reduces to the following integrals after some 
straightforward case distinctions: 

ff dxdy = f 1 ^R 2 -x 2 -y Q dx (C.l) 

J J A Jxq 

= - (x V/? 2 - x 2 + fl 2 arctan ( *' Vl - xy Q 
2 1 X-IW^x 2 )) xn 



